Preliminaries

<<ProPac`

?ProPac`GeoTools`*

Auxillary Functions

Replacement Rules

ReplacementRule[x_List]:=Map[(#->ToExpression[ToString[Head[#]]<>ToString[#[[1]]]])&,x];

ReplacementRule[{V[1],V[2],theta[1]}]

InputReplacementRule[u_List]:=Map[(#->ToExpression[ToString[Head[#]]<>ToString[#[[1]]]])&,u];

InputReplacementRule[{P[1],Q[1]}]

ReplacementRule[{P[1],Q[1]}]

InitialValueRule[x_List]:=Map[(#->ToExpression[ToString[Head[#]]<>"0"<>"["<>ToString[#[[1]]]<>"]"])&,x];

InitialValueRule[{V[1],V[2],theta[1]}]

Newton Raphson

NewtonRaphson[EquilibriumlEquations_List,Variables_List,Parameters_List,InitialVariables_List,InitialParameters_List, ErrorNorm_,MaxIterations_, opts___ ]:= Module[{YY,XX,EE,FF,Δ,InitialInputRules,NewtEqns,sols,k=MaxIterations,i},
InitialInputRules=Inner[Rule,Parameters,InitialParameters,List];
sols=Inner[Rule,Variables,InitialVariables,List];
Δ=Variables;
EE=EquilibriumlEquations/.InitialInputRules;
XX=Chop[Jacob[EE,Variables]];
For[i=1,i≤k,i++,
YY=XX/.sols;
FF=EE/.sols;
NewtEqns=Inner[Equal,YY.Δ-YY.(Δ/.sols)+FF,ConstantArray[0,Length[Δ]],List];
sols=Solve[NewtEqns,Δ][[1]]
];
sols
]

NewtonRaphsonSeydel[EquilibriumEquations_List,Variables_List,Parameters_List,InitialVariables_List,InitialParameters_List,v0_List,λ0_]:= Module[{YY,XX,EE1,EE2,EE3,EE,FF,Δ,Δ0,InitialInputRules,NewtEqns,InitialVariablesRules,InitialParametersRules,V,sols,Data},
InitialVariablesRules=Inner[Rule,Variables,InitialVariables,List];
InitialParametersRules=Inner[Rule,Parameters,InitialParameters,List];
V=Table[v[i],{i,Length[Variables]}];
Δ=Join[Variables,V,Parameters];
Δ0=Join[InitialVariables,v0,InitialParameters];
sols=Inner[Rule,Δ,Δ0,List];
Data={Join[Variables,Parameters]/.sols};
EE1=EquilibriumEquations;
XX=Chop[Jacob[EE1,Variables]];
EE2=(XX-λ*IdentityMatrix[Length[V]]).V;
EE3={V.V-1};
EE=Join[EE1,EE2,EE3];
For[i=1,i≤10,i++,
μ=λ0-i*λ0/5;
sols=NewtonRaphson[EE,Δ,{λ},Δ/.sols,{μ}, 10^-6,5];
Data=Join[Data,{Join[Variables,Parameters]/.sols}]
];
Data
]

Example 1 - 3 Bus, B=0

Set Up Equations

Eqns={-V Sin[θ-ϕ] -V Sin[ϕ]+ΔP==0,V (Sin[ϕ]+Sin[ϕ-θ])-P3==0,-V (Cos[ϕ]+Cos[ϕ-θ]) + (2-B) V^2 - Q3==0};

F={-V Sin[θ-ϕ] -V Sin[ϕ]+ΔP,V (Sin[ϕ]+Sin[ϕ-θ])-P3,-V (Cos[ϕ]+Cos[ϕ-θ]) + (2-B) V^2 - Q3};

FF=F/.{ΔP->0,B->0,Q3->0}

J=Jacob[FF,{θ,ϕ,V}];
J//MatrixForm

NewtonRaphson[FF,{θ,ϕ,V},{P3},{0.0001,-π/4,1/1.38},{-1}, 10^-6,10]

J/.{θ->0,ϕ->-π/4,V->1/(Sqrt[2])}

Eigenvalues[J/.{θ->0,ϕ->-π/4,V->1/(Sqrt[2])}]

Det[J/.{θ->0,ϕ->-π/4,V->1/(Sqrt[2])}]

NullSpace[J/.{θ->0,ϕ->-π/4,V->1/(Sqrt[2])}]

Set Starting Point

Sols=NewtonRaphson[FF,{θ,ϕ,V},{P3},{0.000,0,1},{-0.000}, 10^-6,10]
Data={{θ,ϕ,V,0}}/.Sols
P=0;

NR Until Near Singular Solution

For[j=1,j≤10,j++,
P=-0.1*j;
Sols=NewtonRaphson[FF,{θ,ϕ,V},{P3},{θ,ϕ,V}/.Sols,{P}, 10^-6,10];
Data=Join[Data,{{θ,ϕ,V,P}}/.Sols]
];
Data

Set Up NRS Starting Point

ES=Eigensystem[J/.Inner[Rule,{θ,ϕ,V,P},Data[[-2]],List]]

λ0=ES[[1,-1]]
v0=ES[[2,-1]]

Run NRS Around Nose

Data[[-1]][[{1,2,3}]]
Data[[-1]][[{4}]]

Data1=NewtonRaphsonSeydel[FF,{θ,ϕ,V},{P3},Data[[-2]][[{1,2,3}]],Data[[-2]][[{4}]],v0,λ0]

Continue NR

ResultRules=Map[Inner[Rule,{θ,ϕ,V,P3},#,List]&,Join[Drop[Data,-1],Data1]]

gr1=ListPlot[{-P3,V}/.ResultRules,Joined->True,PlotRange->{{-0.,1.2},{0.4,1.2}},AxesLabel->{"-(Subscript[P, 3])","(Subscript[V, 3])"},AspectRatio->1/2,PlotStyle->Black]

Example 1 - 3 Bus, Q=0.05

Set Up Equations

Eqns={-V Sin[θ-ϕ] -V Sin[ϕ]+ΔP==0,V (Sin[ϕ]+Sin[ϕ-θ])-P3==0,-V (Cos[ϕ]+Cos[ϕ-θ]) + (2-B) V^2 - Q3==0};

F={-V Sin[θ-ϕ] -V Sin[ϕ]+ΔP,V (Sin[ϕ]+Sin[ϕ-θ])-P3,-V (Cos[ϕ]+Cos[ϕ-θ]) + (2-B) V^2 - Q3};

FF=F/.{ΔP->0,B->0.00,Q3->0.05}

J=Jacob[FF,{θ,ϕ,V}];
J//MatrixForm

NewtonRaphson[FF,{θ,ϕ,V},{P3},{0.0001,-π/4,1/1.38},{-1}, 10^-6,10]

J/.{θ->0,ϕ->-π/4,V->1/(Sqrt[2])}

Eigenvalues[J/.{θ->0,ϕ->-π/4,V->1/(Sqrt[2])}]

Det[J/.{θ->0,ϕ->-π/4,V->1/(Sqrt[2])}]

NullSpace[J/.{θ->0,ϕ->-π/4,V->1/(Sqrt[2])}]

Set Starting Point

Sols=NewtonRaphson[FF,{θ,ϕ,V},{P3},{0.000,0,1},{-0.000}, 10^-6,10]
Data={{θ,ϕ,V,0}}/.Sols
P=0;

NR Until Near Singular Solution

For[j=1,j≤10,j++,
P=-0.1*j;
Sols=NewtonRaphson[FF,{θ,ϕ,V},{P3},{θ,ϕ,V}/.Sols,{P}, 10^-6,10];
Data=Join[Data,{{θ,ϕ,V,P}}/.Sols]
];
Data

Set Up NRS Starting Point

ES=Eigensystem[J/.Inner[Rule,{θ,ϕ,V,P},Data[[-2]],List]]

λ0=ES[[1,-1]]
v0=ES[[2,-1]]

Run NRS Around Nose

Data[[-1]][[{1,2,3}]]
Data[[-1]][[{4}]]

Data1=NewtonRaphsonSeydel[FF,{θ,ϕ,V},{P3},Data[[-2]][[{1,2,3}]],Data[[-2]][[{4}]],v0,λ0]

Continue NR

Data1

ResultRules=Map[Inner[Rule,{θ,ϕ,V,P3},#,List]&,Join[Drop[Data,-1],Data1]]

gr2=ListPlot[{-P3,V}/.ResultRules,Joined->True,PlotRange->{{-0.,1.2},{0.4,1.2}},AxesLabel->{"-(Subscript[P, 3])","(Subscript[V, 3])"},AspectRatio->1/2,PlotStyle->{Dashed,Black}]

Show[gr1,gr2]

Example 1 - 3 Bus, Q=-0.05

Set Up Equations

Eqns={-V Sin[θ-ϕ] -V Sin[ϕ]+ΔP==0,V (Sin[ϕ]+Sin[ϕ-θ])-P3==0,-V (Cos[ϕ]+Cos[ϕ-θ]) + (2-B) V^2 - Q3==0};

F={-V Sin[θ-ϕ] -V Sin[ϕ]+ΔP,V (Sin[ϕ]+Sin[ϕ-θ])-P3,-V (Cos[ϕ]+Cos[ϕ-θ]) + (2-B) V^2 - Q3};

FF=F/.{ΔP->0,B->0.00,Q3->-0.05}

J=Jacob[FF,{θ,ϕ,V}];
J//MatrixForm

NewtonRaphson[FF,{θ,ϕ,V},{P3},{0.0001,-π/4,1/1.38},{-1}, 10^-6,10]

J/.{θ->0,ϕ->-π/4,V->1/(Sqrt[2])}

Eigenvalues[J/.{θ->0,ϕ->-π/4,V->1/(Sqrt[2])}]

Det[J/.{θ->0,ϕ->-π/4,V->1/(Sqrt[2])}]

NullSpace[J/.{θ->0,ϕ->-π/4,V->1/(Sqrt[2])}]

Set Starting Point

Sols=NewtonRaphson[FF,{θ,ϕ,V},{P3},{0.000,0,1},{-0.000}, 10^-6,10]
Data={{θ,ϕ,V,0}}/.Sols
P=0;

NR Until Near Singular Solution

For[j=1,j≤10,j++,
P=-0.1*j;
Sols=NewtonRaphson[FF,{θ,ϕ,V},{P3},{θ,ϕ,V}/.Sols,{P}, 10^-6,10];
Data=Join[Data,{{θ,ϕ,V,P}}/.Sols]
];
Data

Set Up NRS Starting Point

ES=Eigensystem[J/.Inner[Rule,{θ,ϕ,V,P},Data[[-2]],List]]

λ0=ES[[1,-1]]
v0=ES[[2,-1]]

Run NRS Around Nose

Data[[-1]][[{1,2,3}]]
Data[[-1]][[{4}]]

Data1=NewtonRaphsonSeydel[FF,{θ,ϕ,V},{P3},Data[[-2]][[{1,2,3}]],Data[[-2]][[{4}]],v0,λ0]

Continue NR

Data1

ResultRules=Map[Inner[Rule,{θ,ϕ,V,P3},#,List]&,Join[Drop[Data,-1],Data1]]

gr3=ListPlot[{-P3,V}/.ResultRules,Joined->True,PlotRange->{{-0.,1.2},{0.4,1.2}},AxesLabel->{"-(Subscript[P, 3])","(Subscript[V, 3])"},AspectRatio->1/2,PlotStyle->{Dotted,Black}]

Show[gr1,gr2,gr3]

Example 4

f={2 Sin[θ1]+Sin[θ1-θ2]+Sin[θ2]-ΔP1,2 Sin[θ2]+Sin[θ2-θ1]+Sin[θ1]-ΔP1}

J=Jacob[f,{θ1,θ2}];

J//MatrixForm

J0=J/.{θ1->π/2,θ2->π,ΔP1->1,ΔP2->2}

NullSpace[{{0,1}}]

NullSpace[{{-1,-2}}]

SpanD[Transpose[J0]]

LeftInverse[J0]

RightInverse[J0]

ff=FullSimplify[{-2,1}.f]

α0=ff/.{θ1->0